Setup

source("helpers.R")

dir.create("results/TLR4/DEA", showWarnings = F, recursive = T)

knitr::opts_chunk$set(fig.width = 10, dpi = 300, results = "hold", fig.show = "hold")
# Heatmaps

hm_cluster_rows <- TRUE # Genes
hm_cluster_cols <- TRUE # Samples
hm_scale_by_row <- TRUE
hm_max_rows <- 30

heatmap.colors <- colorRamp2(c(-2, 0, 2), c("darkblue", "white", "darkred"))

# Volcano plot

vp_lfc_limit <- 5

Supernatant

data <- read_DIA_report("data/TLR4/20240828_120228_ASC_TERT_TLR_10_4_SN_Report.tsv")

data <- data[, data$Cell.Type != "ASC_TERT_TLR10_LOV"]

data_MS1.df <- assay(data, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
data_MS2.df <- assay(data, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")  
data.log2 <- data
assays(data.log2) %<>% lapply(log2)

data_MS1.log2.df <- assay(data.log2, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.log2.df <- assay(data.log2, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")


plot_intensity_boxplot(data_MS1.log2.df, data_MS2.log2.df)

meanSdPlot(assay(data.log2, "MS1"))
meanSdPlot(assay(data.log2, "MS2"))

data.norm <- data
assays(data.norm) %<>% lapply(normalize_matrix, "vsn")

data_MS1.norm.df <- assay(data.norm, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.norm.df <- assay(data.norm, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.norm.df, data_MS2.norm.df)

meanSdPlot(assay(data.norm, "MS1"))
meanSdPlot(assay(data.norm, "MS2"))

plot_pca(data.norm, scale = T, plot_all = T, maxPC = 4)

contrast_list <- c(
  LightVsDark = "ASC_TERT_TLR4_LOV.Light - ASC_TERT_TLR4_LOV.Dark",
  LPSVsDark = "ASC_TERT_TLR4_LOV.LPS - ASC_TERT_TLR4_LOV.Dark",
  LightVsLPS = "ASC_TERT_TLR4_LOV.Light - ASC_TERT_TLR4_LOV.LPS"
)

fit <- fit_DEqMS_model(data.norm, contrast_list)

VarianceBoxplot(fit)

ASC TERT TLR4LOV Light vs. Dark

current_contrast <- 1
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.norm[rownames(data.norm) %in% res.sig$gene, 
                         str_detect(colnames(data.norm), "Light|Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  
print(paste("Found", nrow(data.sig.df), "differentially expressed proteins."))

write.csv(res, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_Light_vs_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_Light_vs_Dark_SN_filtered.csv",
          row.names = FALSE)
## [1] "Found 211 differentially expressed proteins."
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_TERT_TLR4_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TERT TLR4LOV Light vs. Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 180 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TERT TLR4LOV LPS vs. Dark

current_contrast <- 2
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.norm[rownames(data.norm) %in% res.sig$gene, 
                         str_detect(colnames(data.norm), "LPS|Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  
print(paste("Found", nrow(data.sig.df), "differentially expressed proteins."))

write.csv(res, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_LPS_vs_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_LPS_vs_Dark_SN_filtered.csv",
          row.names = FALSE)
## [1] "Found 615 differentially expressed proteins."
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_TERT_TLR4_LOV" = "lightgreen"),
                                            Condition=c("LPS" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TERT TLR4LOV LPS vs. Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 596 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TERT TLR4LOV Light vs. LPS

current_contrast <- 3
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.norm[rownames(data.norm) %in% res.sig$gene, 
                         str_detect(colnames(data.norm), "Light|LPS")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  
print(paste("Found", nrow(data.sig.df), "differentially expressed proteins."))

write.csv(res, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_Light_vs_LPS_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_Light_vs_LPS_SN_filtered.csv",
          row.names = FALSE)
## [1] "Found 597 differentially expressed proteins."
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_TERT_TLR4_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "LPS" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TERT TLR4LOV Light vs. LPS Supernatant", vp_lfc_limit)
## Warning: ggrepel: 565 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

Whole Cell Lysate

data <- read_DIA_report("data/TLR4/20240828_120420_ASC_TERT_TLR_10_4_WCL_Report.tsv")

data <- data[, data$Cell.Type != "ASC_TERT_TLR10_LOV"]

data_MS1.df <- assay(data, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
data_MS2.df <- assay(data, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")  
data.log2 <- data
assays(data.log2) %<>% lapply(log2)

data_MS1.log2.df <- assay(data.log2, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.log2.df <- assay(data.log2, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")


plot_intensity_boxplot(data_MS1.log2.df, data_MS2.log2.df)

meanSdPlot(assay(data.log2, "MS1"))
meanSdPlot(assay(data.log2, "MS2"))

data.norm <- data
assays(data.norm) %<>% lapply(normalize_matrix, "vsn")

data_MS1.norm.df <- assay(data.norm, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.norm.df <- assay(data.norm, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.norm.df, data_MS2.norm.df)

meanSdPlot(assay(data.norm, "MS1"))
meanSdPlot(assay(data.norm, "MS2"))

plot_pca(data.norm, scale = T, plot_all = T, maxPC = 4)

contrast_list <- c(
  LightVsDark = "ASC_TERT_TLR4_LOV.Light - ASC_TERT_TLR4_LOV.Dark",
  LPSVsDark = "ASC_TERT_TLR4_LOV.LPS - ASC_TERT_TLR4_LOV.Dark",
  LightVsLPS = "ASC_TERT_TLR4_LOV.Light - ASC_TERT_TLR4_LOV.LPS"
)

fit <- fit_DEqMS_model(data.norm, contrast_list)

VarianceBoxplot(fit)

ASC TERT TLR4LOV Light vs. Dark

current_contrast <- 1
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.norm[rownames(data.norm) %in% res.sig$gene, 
                         str_detect(colnames(data.norm), "Light|Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  
print(paste("Found", nrow(data.sig.df), "differentially expressed proteins."))

write.csv(res, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_Light_vs_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_Light_vs_Dark_WCL_filtered.csv",
          row.names = FALSE)
## [1] "Found 252 differentially expressed proteins."
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_TERT_TLR4_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TERT TLR4LOV Light vs. Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 228 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TERT TLR4LOV LPS vs. Dark

current_contrast <- 2
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.norm[rownames(data.norm) %in% res.sig$gene, 
                         str_detect(colnames(data.norm), "LPS|Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  
print(paste("Found", nrow(data.sig.df), "differentially expressed proteins."))

write.csv(res, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_LPS_vs_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_LPS_vs_Dark_WCL_filtered.csv",
          row.names = FALSE)
## [1] "Found 318 differentially expressed proteins."
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_TERT_TLR4_LOV" = "lightgreen"),
                                            Condition=c("LPS" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TERT TLR4LOV LPS vs. Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 286 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TERT TLR4LOV Light vs. LPS

current_contrast <- 3
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.norm[rownames(data.norm) %in% res.sig$gene, 
                         str_detect(colnames(data.norm), "Light|LPS")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  
print(paste("Found", nrow(data.sig.df), "differentially expressed proteins."))

write.csv(res, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_Light_vs_LPS_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/TLR4/DEA/ASC_TERT_TLR4LOV_Light_vs_LPS_WCL_filtered.csv",
          row.names = FALSE)
## [1] "Found 177 differentially expressed proteins."
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_TERT_TLR4_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "LPS" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TERT TLR4LOV Light vs. LPS Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 152 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps